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ABSTRACT 


Parametric  piecewise  cubic  polynomials  are  used  throughout  the  computer 
graphics  industry  to  represent  geometric  curved  shapes.  The  exploration  of  the  use  of 
parametric  curves  and  surfaces  can  be  viewed  as  the  birth  of  Computer  Aided 
Geometric  Design  (CAGD).  In  this  thesis,  least  squares  approximation  is  used  for 
fltting  a  geometrically  continuous  (G')  piecewise  parametric  cubic  polynomial  to  a 
sequence  of  ordered  points  in  the  plane.  Cubic  Bdzier  curves  are  used  as  a  basis. 

The  parameterization,  the  control  points,  the  number  of  knots,  and  their  locations  are 
determined  as  part  of  the  approximation  process.  A  development  of  the  algorithm  is 
given,  along  with  some  results  for  a  variety  of  sets  of  ordered  data. 
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I.  INTRODUCTION 


A.  GENERAL 

In  many  CAGD  applications,  a  user  wishes  to  produce  a  smooth  curve  from  a 
given  ordered  set  of  data  points  such  that  the  results  captures  the  general  shape  of  the 
data.  This  is  no  easy  problem  as  many  numerical  analysts  have  discovered.  There 
are  several  types  of  approximations  that  can  be  used  to  solve  this  problem.  In  simple 
cases  polynomial  interpolation  is  generally  used.  In  advanced  cases,  piecewise 
polynomials  are  used.  These  methods  are  incorporated  into  least-squares  data  fitting 
in  approximating  curved  shapes  from  a  set  of  discrete  data  points.  When  parametric 
curves  are  used,  an  appropriate  parameterization  is  essential  in  obtaining  a  good 
representation  of  the  original  shape. 

The  method  of  fitting  given  data  by  a  polynomial  which  coincides  with  the 
data  at  certain  specified  points  is  in  some  cases  a  rather  inefficient  one.  In  particular, 
when  the  function  is  known  at  all  points  in  an  interval,  it  seems  undesirable  to  select 
only  a  relatively  small  set  of  arbitrarily  chosen  points  at  which  to  "match",  Weeg 
[Ref.  1:  p.  126].  Obviously,  using  a  large  number  of  points  results  in  a  prohibitively 
high  degree  polynomial  for  numerical  computations.  Another  undesirable  case  is  one 
in  which  the  degree  of  reliability  of  such  values  is  not  clearly  established,  as  is  often 
the  case  for  experimental  data. 

The  method  of  least  squares  is  designed  to  treat  these  two  classes  of  problems. 
The  usual  criterion  of  least-squares  approximation  is  to  minimize  the  sum  of  the 
squares  of  the  errors  within  a  specified  domain.  Detailed  explanations  of  least 
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squares  data-fitting  can  be  found  in  most  books  on  "Numerical  Analysis",  [Ref.  l-S], 

Polynomial  interpolation  is  the  most  fundamental  of  all  approximation 
concepts.  The  algebraic  polynomials  p,(x)  are  by  far  the  most  important  and  popular 
approximating  functions,  Carnahan  [Ref.  2:  p.  2].  The  theory  of  polynomial 
interpolation  is  well  developed  and  fairly  simple.  Polynomials  are  easy  to  evaluate 
and  their  sums,  products,  and  differences  are  also  polynomials.  Polynomials  can  be 
differentiated  and  integrated  with  little  difficulty.  In  addition,  if  the  origin  of  the 
coordinate  system  is  shifted  or  if  the  scale  of  the  independent  variable  is  changed,  the 
transformed  polynomials  remain  polynomials,  that  is,  if  p,(x)  is  a  polynomial,  so  are 
pa(x+a)  and  PB(ax),  Ralston  [Ref.  3:  p.  34]. 

All  these  advantages  of  the  polynomial  would  be  of  little  value  if  there  were 
no  analytical  justification  for  believing  that  polynomials  can,  in  fact,  yield  good 
(implies  the  error  can  be  made  arbitrarily  small)  approximations  for  a  given  function 
f(x).  This  theoretical  justification  does  exist;  any  continuous  function  f(x)  can  be 
approximated  to  any  desired  degree  of  accuracy  on  a  specified  closed  interval  by 
some  polynomial  p„(x).  This  follows  from  the  Weierstrass  approximation  theorem 
stated  here: 

If  f(x)  is  continuous  in  the  closed  interval  [a,b],  then  given  any  6>0,  there  is 
some  polynomial  Pa(x)  of  degree  n3n(c)  such  that 

|f(x)-p.(x)|<6,  X€[a,b]. 

It  is  important  to  note  that  the  usual  criteria  for  generating  interpolating  polynomials 
in  no  way  guarantees  that  the  polynomial  found  is  the  one  which  the  Weierstrass 
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theorem  shows  must  exist.  This  is  why  polynomial  interpolation  is  mostly  of 
theoretical  value,  while  faster  and  more  accurate  methods,  piecewise  polynomials,  are 
used  extensively  in  approximation  analysis. 

Piecewise  polynomials  are  more  widely  used  in  least-squares  data  fitting  since 
they  are  more  flexible  than  simple  polynomials.  If  we  wish  to  approximate  to  a 
function  f  on  [a,b],  we  may  divide  [a,b]  into  N  sub-intervals,  by  introducing  points  of 
sub-division 

a  =  Xo<x,<  .  .  .  <  x„  =  b, 

and  use  some  form  of  approximating  or  interpolating  polynomial  on  each  sub-interval 
[Xi.i,xJ,  I  ^i^N.  This  is  called  piecewise  polynomial  approximation.  Phillips  [Ref. 
4:  p.  113]. 

Generally,  the  approximating  polynomials  used  on  neighboring  intervals 
[Xi.i,Xi]  and  [xj,  Xj+J  will  not  take  the  same  value  at  the  common  point  Xj.  It  is 
natural  to  seek  piecewise  approximations  which  are  continuous  and  whose  first  few 
derivatives  are  also  continuous  at  the  nodes  of  subdivision,  so  that  the  connection 
appears  smooth.  These  are  called  spline  approximations.  The  most  common 
piecewise  polynomial  approximation  using  cubic  polynomials  between  each  successive 
pair  of  nodes  is  cubic  spline  interpolation.  A  piecewise  cubic  polynomial 
approximation  to  a  curved  shape  consists  of  a  number  of  single  cubic  segments 
connected  end-to-end  to  form  a  continuous  curve.  A  general  cubic  polynomial 
involves  four  constants;  so  there  is  sufficient  flexibility  in  the  cubic  spline  procedure 
to  ensure  not  only  that  the  interpolant  is  continuous  on  the  interval,  but  also  that  it  has 


3 


a  continuous  first  and  second  derivative  on  the  interval.  But,  as  stated  by  Burden 
[Ref.  5:  p.  132]  the  constmction  of  the  cubic  spline  does  not  assume  that  the 
derivatives  of  the  in^erpolant  agree  with  those  of  the  function,  even  at  the  nodes. 

If  X  and  y  are  given  as  functions 

X  =  f(t)  y  =  g(t) 

over  an  interval  of  t-values,  then  the  set  of  points  (x,y)  ==  (f(t),g(t))  defined  by  these 
equations  is  called  a  curve  in  the  coordinate  plane.  The  equations  are  parametric 
equations  for  the  curve.  The  variable  t  is  called  the  parameter  of  the  curve.  When 
we  write  parametric  equations  for  the  curve  in  the  plane,  we  say  that  we  have 
parameterized  the  curve.  For  example,  a  parametric  equation  of  the  unit  circle  is  as 
follows: 

X  =  cos  t  y  =  sin  t  0:Sts2T  . 

The  spline  interpolation  problem  is  usually  stated  as  "given  data  points  q^  and 
parameter  values  Uj,...".  The  simplest  way  to  determine  the  Uj  is  to  set  Uj  =  i.  This 
is  called  uniform  or  equidistant  parameterization,  Farin  [Ref.  6:  p.  130].  This 
method  is  not  useful  in  most  practical  situations.  The  reason  is  that  uniform 
parameterization  does  not  take  into  account  the  geometry  of  the  data  points,  for  an 
heuristic  explanation  see  Farin  [Ref.  6:  p.  130].  There  are  exceptions  in  which 
uniform  parameterization  fares  better  than  other  methods,  according  to  Foley  [Ref.  7: 
p.  86]. 

The  most  commonly  used  parameterization  or  knot  spacing  is  chord 
length  spacing,  where  hj=  Iq^^,  -  q^j  (Euclidean  distance).  Some  authors  have 
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suggested  using  chord  length  because  it  approximates  the  arc  length  of  a  parametric 
curve,  for  example  Ahlberg  [Ref.  8:  pp.  2S4-2S8].  Another  advantage  for  chord 
length  parameterization,  as  proved  by  Epstein  [Ref.  9:  pp.  261-268],  is  that  it 
guarantees  there  will  be  no  cusps  for  the  case  of  a  closed  periodic  curve.  When  the 
data  is  poorly  scaled  or  when  the  direction  of  the  data  changes  abruptly,  the  chord 
length  knot  spacing  often  produces  visually  poor  results,  as  seen  in  Foley  and  Nielson 
[Ref.  10:  pp.  261-271].  A  curve  may  be  parametrized  in  many  ways  and  the  choice 
of  parametization  is  important  as  discussed  by  de  Boor  [Ref.  11:  pp.  315-31 8].  There 
is  probably  no  "besf*  parametrization  since  most  known  methods  can  be  defeated  by  a 
suitably  chosen  data  set.  For  our  algorithm,  chord  length  parametrization  was  chosen 
as  an  initial  guess. 

The  problem  presented  here  is  the  following: 

Given  an  ordered  set  m  of  data  points  Q  =  (Xj,yi)  for  i  =  l,...,m,  we  wish  to  fit  a 
piecewise  Bdzier  cubic  polynomial  B(t)  =  (x(t),y(t))  to  Qi  that  minimizes  the  sum  of 
squares  of  distances  from  the  data  points  to  the  nearest  point  on  the  piecewise  B6zier 
curve.  The  optimization  routine  "fmins"  minimizes  this  sum  by  optimizing  the  free 
parameters  (i.e.,  knot  points,  control  points,  unit  tangent  vectors).  In  particular  we 
wish  to  fit  geometrically  continuous  (G')  curves,  rather  than  parametrically  continuous 
(C‘)  curves,  generally  because  the  former  is  a  less  restrictive  condition. 


B.  CONTINUITY/PARAMETRIC  AND  GEOMETRIC 


Parametric  continuity  (CO  of  a  curve  means  that  the  component  functions  are  r 
times  continuously  differentiable  in  a  given  interval  [a,b].  Therefore,  C  continuity 
implies  that  the  component  functions  are  one  time  continuously  differentiable. 
Furthermore,  C  continuity  is  based  on  the  interplay  between  domain  and  range 
conHgurations.  Without  any  information  on  the  domain  of  a  curve,  we  cannot  make 
any  statements  concerning  differentiability.  In  another  case,  zero  tangent  vectors  may 
give  rise  to  comers  or  cusps  in  curves  (which  are  C),  a  fact  that  intuitively 
contradicts  the  concept  of  differentiability.  Smoothness  and  differentiability  only 
agree  for  functional  curves,  the  connection  between  them  is  lost  in  the  parametric 
case.  Differentiable  curves  may  not  be  smooth  and  smooth  curves  may  not  be 
differentiable.  We  emphasize  that  by  "smooth",  we  mean  "perceptually  smooth". 

Geometric  continuity  of  order  r  for  a  curve  implies  that  the  curve  is  r  times 
differentiable  with  respect  to  arc  length,  Farin  [Ref.  6].  It  also  concerns  how 
parametric  curves  and  surfaces  can  be  pieced  together  in  a  smooth  way.  If  a  curve 
has  continuous  curvature,  it  is  (second  order  geometrically  continuous).  Curves 
with  a  zero  tangent  vector  under  the  given  parametrization,  can  be  G^  for  example 
the  parametric  curve  x  =  t*,  y  =  t’  is  G^  at  t  =  0.  A  tutorial  view  of  the  concepts  of 
geometric  continuity  within  the  subject  of  CAGD  is  presented  by  Gregory  [Ref.  12: 
pp.  353-371].  The  following  examples  illustrate  the  difference  between  these  two 
concepts; 
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(1)  Consider  the  parametric  curve  defined  by 


x  =  e.  y=|t>|, 

^  .  dim  .  3 1 1 1  ‘  .  3 1 1 M '-rlT  •  3 1 1 1 1: . 

dt  dt  '  '  dt  '  '  Th 

Note  that  y  =  |x| .  This  curve  is  C*  but  not  G‘,  see  Figure  1. 

(2)  Consider  the  Bdzier  parametric  curves  for  the  control  points  (0,0),  (0,3),  (3,3), 
and  (3,3),  (6,3),  (6,0)  respectively. 


(a) 

X 

II 

y(t)  =  6t-3t^ 

t  6  [0,1];  and 

(b) 

x(t)  =  12t-3t^-6 

y(t)  =  6t-3t^ 

[1,2]. 

For  curve  (a). 


-^=6t 

dt 


^=6-6t, 

dt 


For  curve  (b). 


and  is  same  as  (a)  . 
dt  dt^ 


Thus,  y(t)  is  infinitely  differentiable,  but  x(t)  is  not  twice  differentiable  at  t=l;  as 
shown  above  x^fl)  =  6  for  (a)  and  x"  (1)  =  -6  for  (b).  However,  the  curve  does 
have  continuous  curvature  at  t  =  1,  as  shown  below.  The  curvature  of  a  plane  curve  is 
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ngure  1.  X =1^,  y  =  1 I .  This  curve  is  C,  thus  parametrically  continuous  but  not 
G‘. 
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defined  as  the  following: 


/  dx\i  dM  _  f  dy\(  d^x 

m  *  mf 

Therefore,  for  (a), 

rr  ^  (6)  (-6)  -  (6-6)  (6)  ^  -36  _  -1 

1  216  6  ' 
[(6)2  +  (6-6)2]  2 


rr  ^  (12-6)  (-6)  -  (6-6)  (6)  ^  -36  _  -1 

1  216  6  ' 
[(12-6)2  +  (6-6)2]  2 

Both  curves  have  curvature  equal  to  *1/6.  Figure  2  depicts  a  graph  of  these  two 
curves  that  coincide  at  point  (3,3). 
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0  1  2  3  4  5 

X 

Figure  2.  This  curve  consists  of  two  piecewise  parabolic  arcs,  that  coincide  at  P2. 
This  curve  is  not  twice  differentiable  at  P2,  but  it  has  continuous  curvature,  and  is 
therefore  G^. 
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C.  BERNSTEIN  POLYNOMIALS 


Bernstein  polynomials  are  a  useful  tool  for  matching  a  curve  to  a  function  f  on 
a  closed  bounded  interval  [a,b].  The  polynomials  are  so  named  because  of  their 
central  role  in  Bernstein’s  constructive  proof  of  the  Weierstrass  approximation 
theorem.  The  essence  of  this  proof  is  the  construction  of  a  sequence  of  polynomials 
which  can  be  shown  to  converge  uniformly  to  f.  It  also  turns  out  that  the  derivatives 
of  these  approximations  converge  uniformly  to  those  of  f  -  assuming  the  latter  exists. 
Buchanan  [Ref.  13]. 

We  begin  with  the  definition  of  the  Bernstein  polynomials. 

Definition:  The  Bernstein  polynomial  of  degree  n  associated  with  the  function 
f  on  [a,b]  is  defined  by 

jTo  '  ' 

where  the  points  Xj  =  a  +  ih  =  a  +  (i/n)(b-a)  for  i  =  0,l,...,n.  In  the  special  case 
where  the  interval  is  [0,1],  equation  (1)  reduces  to 

jZo  '  '  ** 

where  the  Bernstein  basis  polynomials  are  defined  by 

Ba,i(x)  =  (i  =  0, 1,  .  .  .  n  =  0, 1, 2,  .  .  . )  . 

Those  who  have  studied  probability  theory  should  recognize  the  Bernstein 
basis  polynomials  as  the  probability  density  functions  for  a  binomial  distribution. 
Specifically,  Bg,j(x)  is  the  probability  of  achieving  exactly  i  successes  in  a  sequence  of 


n  independent  trials  where  the  probability  of  success  on  any  one  trial  is  x.  Bernstein 
worked  in  probability  theory  as  well  as  in  differential  equations  and  approximation 
theory. 


In  view  of  this,  it  is  clear  that 

=  (x+(l-x))«  =  1  (3) 

from  which  it  can  be  deduced  that  the  Bernstein  polynomial  for  the  function  f(x)  ■  1  is 
itself  identically  1;  equation  (3)  is  the  "partition  of  unity  property”. 

As  another  example,  by  differentiating  (3)  we  obtain 


0  * 

and  multiplying  by  x(l-x)  yields 


-  (n-i)x] 


V  (5)x  ^(l-x)  ( i  -nx) 

fZo  '  ' 


=  0 


Therefore,  it  follows  that 


=  X  (4) 

using  (3).  Thus,  the  Bernstein  polynomial  for  the  function  f(x)=x  reproduces  this 
function  exactly. 

Continuing  in  this  way,  differentiating  (4)  we  have 
and  multiplying  by  x(l-x)/n  yields 
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1 


(1-x) 


(i-nx) 


or 


Since 


x(l-x) 

n 


x^(l-x) 


0-i 


max  x(l  -  x)  =  7 
o<«<i  '  '  4 


it  follows  that 

|x^-S„(x*;x)  L  s  -^-0  as  n-«» 

The  good  j^proximation  properties  of  Bernstein  polynomials  give  us  hope  that  B,(f;x) 
should  be  a  good  approximation  to  f(x)  over  [0, 1].  As  it  turned  out  however, 
Bernstein  polynomials  are  good  at  mimicking  the  behavior  of  f,  but  slow  in 
converging  to  f. 

Bernstein  basis  polynomials  are  usually  used  to  express  Bdzier  curves.  The 
polynomial  is  defined  explicitly  by 
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where  t  c  [0,1]  and  the  binomial  coefficients  are  given  by 


n! 

i !  (ji-i)  ! 


if  Oiiin 


V  0  else 

One  important  property  of  Bernstein  polynomials  is  that  they  satisfy  the  recursion: 


B,i(t)  =  (l-t)B,,,(t)  +  tB,„.,,(t) 

with  Boo(t)B  1  and  B„j(t)*0  for  j  e  {0,1,..., n}.  The  proof  can  be  found  in  Farin 
[Ref.  6].  Also  note  that  Bernstein  polynomials  are  nonnegative  over  the  interval 
[0,1];  this  leads  to  the  convex  hull  property.  The  convex  hull  property  states  that  the 
curve  is  contained  within  the  convex  hull  of  the  control  polygon.  The  convex  hull  of 
a  set  of  points  is  the  smallest  convex  set  containing  the  set  of  points.  Bernstein 
polynomials  are  also  symmetric  with  respect  to  t  and  (1-t),  and  thus  Bgj(t)  = 
Bo,a.j(l-t).  The  Bernstein  polynomial  Bg,j  has  only  one  maximum  and  attains  it  at 
t=i/n;  given  a  function  f(x),  if  the  f(i/n)  value  is  changed,  then  the  Bernstein 
polynomial  is  mostly  affected  by  this  change  in  the  region  of  the  curve  around  the 
parameter  value  t  =  i/n,  which  relates  to  the  local  control  property  exhibited  by 
B-spline  curves. 
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D.  B£zI£R  curves 

Bezier  curves  and  ""rfaces  were  independently  developed  by  P.  de  Casteljau  at 
Citroen  and  by  P.  Bdzier  at  Renault  Automobile  Company.  Bdzier  developed  them  in 
the  early  1960’s  to  fill  a  need  for  curves  whose  shape  can  be  readily  controlled  by 
changing  a  few  parameters  and  which  exhibit  a  "local  control"  property.  Bdzier’s 
application  was  to  construct  pleasing  surfaces  for  car  bodies,  de  Casteljau ’s 
development  was  never  published,  and  so  the  whole  theory  of  parametric  polynomial 
curves  and  surfaces  in  Bernstein  form  now  bears  Bdzier’s  name.  Farin  [Ref.  6]. 

Bdzier  curves  are  numerically  the  most  stable  of  all  polynomial  bases  currently 
used  in  Computer  Aided  Design  (CAD)  systems.  Therefore,  Bdzier  curves  are  the 
ideal  geometric  standard  for  representation  of  piecewise  polynomial  curves.  Bdzier 
curves  also  lend  themselves  easily  to  a  geometric  understanding  of  many  CAGD 
phenomena  and  may  be  used  to  derive  the  theory  of  rational  and  non-rational  B-spline 
curves. 

Curve  fitting  has  received  a  fair  amount  of  attention,  even  before  computer 
graphics  came  along.  Piecewise  polynomial  curves  have  been  used  to  mathematically 
interpolate  discrete  data  sets.  The  algorithm  developed  fits  a  piecewise  cubic  Bdzier 
curve  to  a  set  of  data  points.  In  particular  it  fits  geometrically  continuous  (C) 
curves,  rather  than  parametrically  continuous  (C)  curves. 

The  control  points  or  Bdzier  points  are  coordinate  pairs  (Xj,  yj  for  i=0,l,...N 
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which  are  used  to  define  a  Bdzier  curve  by  setting 


S  H 

x(c)  y(t)  . 

The  polygon  formed  by  (Xo,yo),  (x„  yi),...(x„yj  is  called  the  Bdzier  polygon  or 
control  polygon. 

Suppose  we  are  given  a  set  of  data  points,  p;,  and  the  points  do  not  necessarily 
progress  from  left  to  right.  The  coordinates  of  each  point  are  treated  as  a  two  vector 
components. 


In  parametric  form,  the  curve  is 

The  n-th  degree  Bdzier  polynomial  determined  by  n+1  points  is  given  by 

For  n=2,  this  would  give  the  quadratic  curve  defined  by  three  points  pb.P1.P2: 

P(t)  =  (l)(l-t)^o  +  2(l-t)(t)p,  +  (l)t^  . 

The  parametric  form  of  the  above  curve  is 

x(t)=(l-t)^Xo  +  2(l-t)(t)x,  +  Fxj  , 
y(t)=(l-t)^yo  +  2(l-t)(t)y,  +  Fyz  • 

Note  that  if  t=0,  x(0)=x<,  and  y(0)=yo.  If  t=l,  the  point  is  (X2,y2).  As  t 
takes  on  values  between  0  and  1,  a  curve  is  traced  that  goes  from  the  first  point  to  the 
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third  point  of  the  set.  Usually,  the  curve  will  not  pass  through  the  central  point  of  the 
three  unless  the  points  are  collinear;  if  the  points  are  collinear,  the  curve  is  a  straight 
line  through  all  three  points.  In  effect,  the  points  of  the  second-degree  Bdzier  curve 
have  coordinates  that  are  weighted  sums  of  the  coordinates  of  the  three  points  that  are 
used  to  define  it. 

For  n=3,  we  get  the  cubic  Bdzier  polynomial,  which  we  will  use  in  our 
algorithm.  The  basic  properties  of  higher  degree  Bdzier  polynomials  are  the  same  as 
for  the  cubic.  The  parametric  form  for  the  Bdzier  cubic  is  the  following: 

x(t)=(l-tyxo  +  3(l-t)^(t)x,  +  3(l-t)(F)x2  +  t^xj  , 
y(t)=(l-tyyo  +  3(l-t)2(t)y,  -I-  3(l-t)(e)y2  t^y,  . 

Note  again  that  (x(0),  y(0))  =  po  and  (x(l),y(l))  =  pj,  and  that  the  curve  will  not 
ordinarily  go  through  the  intermediate  points.  As  illustrated  in  Figure  3(a)  and  3(b), 
changing  the  intermediate  "control"  points  changes  the  shape  of  the  curve.  Figure 
3(c)  is  a  sixth  degree  Bdzier  curve,  notice  that  the  curve  tends  to  follow  along  the 
path  of  the  points  (this  will  be  further  explained  later). 

Sometimes  it  is  convenient  to  represent  the  Bezier  curve  in  matrix  form.  For 
Bdzier  cubics,  this  is 


P(t) 


4  t,l] 

O 


-1 

3 

-3 

ijPo 

3 

-6 

3 

o]Pi 

-3 

3 

0 

olpa 

1 

0 

0 

0|p3, 

Some  properties  of  B6zier  cubics  are: 
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(i)  Since 


=  3{Xi-Xo)  and  =  ^(y^-Yo)  at  t*0  , 

the  slope  of  the  curve  at  t=0  is 

dx  (Xj-Xq) 

which  is  the  slope  of  the  secant  line  between  po  and  p,.  Similarly,  the  slope  at  t=l  is 
the  same  as  the  secant  line  between  the  last  two  points. 

(ii)  The  Bdzier  curve  is  contained  in  the  convex  hull  determined  by  the  four 

points. 

If  a  curve  has  a  complex  shape,  then  its  Bdzier  representation  will  have  a  very 
high  degree  (for  practical  purposes,  degrees  exceeding  10  are  prohibitive).  However, 
such  complex  curves  can  be  modeled  using  piecewise  or  composite  Bdzier  curves. 

Bdzier  cubic  curves  can  be  continued  beyond  the  first  set  of  four  points  by 
subdividing  seven  points  (p^  to  p^)  into  two  groups  of  four  points,  with  the  central 
point  Pa  belonging  to  both  sets.  Because  the  curve  is  tangent  to  the  control  polygon  at 
the  end  control  points,  the  piecewise  curve  has  continuous  slope  provided  the  points 
P2,  Pa,  and  P4  are  collinear.  See  Figure  4  for  an  example  of  this. 

A  piecewise  Bdzier  curve  s  is  the  continuous  map  of  a  collection  of  intervals 
Uo<...  <Ul  into  R“,  where  each  interval  [Ui,Ui+,)  is  mapped  onto  a  polynomial  curve 
segment.  Each  real  number  Uj  is  called  a  breakpoint  or  knot.  The  collection  of  all  n, 
is  called  the  knot  sequence.  For  every  parameter  value  u,  there  is  a  corresponding 
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1  2  3  4  5  6  7  8 

X 

Figure  4.  O  control  points,  ®  -*  knot  points,  +  -*  data  points.  The  curve  is 
G‘  continuous  at  P3  because  of  tangency  at  the  end  control  points. 
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point  s(u)  on  the  curve  s.  Let  value  v  be  from  an  interval  [Ui.Ui^J.  We  can  introduce 
a  local  coordinate  (or  local  parameter)  t  for  the  interval  [Ui,Ui+,]  by  setting 


t  =  - —  . 

One  can  check  that  t  varies  from  0  to  1  as  v  varies  from  Uj  to  Uj^.!.  The  collection  of 
the  Bdzier  polygons  for  all  curve  segments  itself  forms  a  polygon;  it  is  called  the 
piecewise  Bdzier  polygon  of  s.  Farin  [Ref.  6]. 

It  is  the  geometry  of  Bdzier  curves  that  makes  them  a  good  choice  for  our 
purposes.  Sometimes,  the  shape  of  the  curve  may  require  a  cusp  or  comer;  this  can 
be  achieved  by  having  knot  points  coincide,  which  produced  zero  tangent  vectors. 
Piecewise  Bdzier  curves  have  the  geometric  property  that  by  changing  one  of  the 
points  we  only  change  one  portion  of  the  curve,  a  "local"  effect.  Generally,  cubic 
spline  curves  have  a  "global"  effect,  in  that  the  curve  from  the  first  to  the  last  point 
would  be  affected.  Bdzier  curves  are  not  really  interpolating  splines,  since  the  curves 
do  not  normally  pass  through  all  of  the  points.  In  this  respect  they  show  some 
similarity  to  least-squares  curves  and  thus  are  effective  in  approximating  data.  Other 
properties  of  Bdzier  curves  are  discussed  by  Farin  [Ref.  6]  and  Gerald  [Ref.  14:  pp. 
217-222]. 


23 


n.  THE  PROBLEM:  HTTING  ORDERED  DATA  IN  A  PLANE 
A.  USING  B£zIER  curves  AS  A  BASIS 

It  is  interesting  to  note  how  piecewise  cubic  polynomials  relate  to  a  more 
common  specification,  Bezier  curves.  In  this  representation,  the  variables  are  vector 
values  called  control  points  which  are  multiplied  by  real-valued  functions  of  t  to  fix 
the  shape  of  the  curve.  The  control  points  are  therefore  coordinate  pairs  (Xj,yi)  for 
i=:0,l,...,N  which  are  used  to  define  a  parametric  curve  -  a  Bdzier  curve  -  by  setting 

n 

xu)  =  ^{t) 

1^0 

y(t)  . 

Such  a  specification  can  be  converted  to  a  scheme  of  vector-valued  basis 
functions  and  real-valued  coefficients  as  follows:  let  the  control  points  be  specified  as 

each  control  point  having  its  own  coordinate  system  with  basis 

Furthermore,  let  the  original  basis  functions  be  {^j(t)}.  Then  the  parametric  curve  is 

=  ^  ( t) )  +  ajj  ( t) ) 
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so  the  sqppropriate  set  of  vector  valued  basis  functions  is  {by^j(t)|i==  1,2}.  Plass  and 
Stone  [Ref.  IS]. 

This  means  we  can  preset  linear  constraints  on  the  control  points  for  the 
resulting  curve.  Since  each  control  point  has  been  given  its  own  coordinate  system, 
by  a^ropriate  choice  of  the  coordinate  systems  and  the  free  coefficients,  each  control 
point  can  be  constrained  to  be  on  a  particular  point  or  line,  or  to  be  free  to  move 
anywhere  in  the  plane. 

It  is  the  control  points  as  seen  in  Figure  1(a)  and  1(b),  that  control  the 
behavior  of  the  Bdzier  curve,  and  in  some  sense  make  the  curve  "mimic”  the  Bdzier 
polygon.  This  is  the  reason  why  Bdzier  curves  provide  such  a  handy  tool  for  the 
design  of  curves:  In  order  to  reproduce  the  shape  of  a  hand  drawn  curve,  it  is 
sufficient  to  specify  a  control  polygon  that  "exaggerates"  the  shape  of  the  curve. 
Letting  the  computer  draw  the  Bdzier  curve  defined  by  the  polygon,  and  if  necessary, 
adjusting  the  location  and  possibly  the  number  of  the  polygon  vertices,  one  will 
reproduce  a  given  curve  after  two  to  three  iterations,  Farin  [Ref.  6]. 
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B.  USING  PIECEWISE  PARAMETRIC  CUBICS  WITH  G‘  CONTINUITY 

Our  design  methodology  is  to  represent  all  shapes  analytically  using  piecewise 
parametric  cubic  polynomials,  from  an  ordered  set  of  discrete  data  points.  A 
parametric  cubic  polynomial  is  quite  flexible,  but  has  some  limits.  For  instance,  it 
can  contain  at  most  one  loop  or,  if  it  has  no  loop,  at  most  two  inflection  points. 
Therefore,  the  least-squares  fit  is  very  vulnerable  to  a  poorly  chosen  continuity 
constraint  such  as  a  bad  tangent  vector. 

Following  Plass  and  Stone  [Ref.  IS],  the  common  specification  of  a  parametric 
cubic  polynomial  is  F(X(t),Y(t)),  where  X  and  Y  are  cubic  polynomials  and  t  lies  in 
the  range  0  to  1.  Using  continuous,  piecewise  functions  will  enable  us  to  constrain 
the  endpoints  and  tangents  for  each  cubic  piece.  To  obtain  G'  continuity  the  knot 
point  must  lie  between  two  adjacent  control  points  on  the  same  line,  as  in  Figure  S;  Pi 
is  a  knot  point.  A  cusp  will  occur  if  the  two  control  points  both  precede  or  follow  the 
knot  point  along  the  same  line.  A  cusp  is  also  defined  as  a  point  on  a  curve  where 
the  first  derivative  vector  vanishes.  Thus  zero  tangent  vectors  may  give  rise  to 
comers  or  cusps  in  curves,  occurring  when  two  or  all  three  points  coincide.  A  comer 
is  a  point  on  a  curve  where  the  tangent,  not  necessarily  the  tangent  vector,  changes  in 
a  discontinuous  way. 
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Figure  5.  O  -*  control  points,  ®  -*  knot  points,  +  data  points.  Shows  curve 
with  G*  continuity  at  the  second  knot  point. 
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C.  ALLOWING  FREE  PARAMETERS 

We  have  seen  the  parametric  form  of  a  Bezier  curve.  B6zier  curves  were 
developed  to  fill  a  need  for  curves  whose  shape  can  be  readily  controlled  by  changing 
a  few  parameters.  By  Figure  6,  one  can  easily  determine  the  parameters  of  the  curve. 
The  parameters  are  the  locations  of  the  knot  points,  the  direction  or  angle  of  the  unit 
tangent  vectors  at  the  knot  points,  and  the  distance  from  the  knot  points  to  adjacent 
control  points.  Note  that  the  unit  tangent  vector  is  parallel  to  the  three  collinear 
control  points.  These  three  sets  of  parameters  determine  the  shape  of  the 
approximating  curve. 

The  number  and  placement  of  the  knot  points  are  critical  in  getting  a  good 
representation  of  the  data,  since  the  approximating  curve  must  pass  through  each  knot 
point.  Choosing  the  knot  locations  can  be  and  has  been  done  in  several  ways.  The 
technique  described  by  Reeves  [Ref.  16]  in  fitting  cubic  pieces  to  a  set  of  data  points 
is  stated  thusly;  "one  simple  method  for  defining  the  knots  is  to  "grow"  the  pieces  out 
until  the  fit  for  that  piece  exceeds  some  threshold.  In  other  words,  starting  with  the 
previous  knot,  keep  adding  data  points  until  the  piece  is  as  long  as  possible.  Each 
new  piece  must  be  constrained  to  maintain  continuity."  This  will  work  in  the  sense 
that  the  resulting  curve  will  fit  everywhere  within  some  specified  tolerance,  but  is 
vulnerable  to  local  phenomena  such  as  a  badly  defined  tangent.  Another  choice,  also 
described  by  Reeves,  is  subdivision. 

Plass  and  Stone  [Ref.  IS]  used  "dynamic  programming"  to  search  a  subset  of 
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Perturbed  Exponential  Spiral 


Hgure  6.  O  -*  control  points,  ®  -*  knot  points,  +  -*  data  points.  Parameters: 
Location  of  the  knots,  the  direction  or  angle  of  the  unit  tangent  vectors,  and  the 
distance  from  the  knot  points  to  adjacent  control  points. 
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the  sample  points  as  potential  knot  positions.  The  criterion  of  Powell  [Ref.  17;  pp. 
65-83]  called  the  "trend"  for  data  fitting,  is  that  knots  are  inserted  successively  until  a 
satisfactory  fit  is  obtained.  An  advantage  of  this  criterion  is  that  it  is  applicable  even 
if  the  magnitude  of  error  is  unknown.  This  method  was  also  used  by  Ichida  and 
Kiyono  [Ref.  18:  pp.  164-174]  in  their  one-pass  method  for  curve  fitting. 

The  unit  tangent  vectors  are  calculated  at  each  knot  point.  The  direction  or 
angle  of  the  unit  tangent  vectors  determined  the  direction  in  which  the  approximating 
curve  will  travel.  They  are  also  used  to  calculate  the  positions  of  the  interior  control 
points  for  each  cubic  segment.  These  control  points  lie  along  the  same  line  as  the 
knot  point  and  their  location  determines  the  geometric  continuity  of  the  curve  at  that 
knot  point;  this  line  is  parallel  to  the  unit  tangent  vector. 

So,  the  collinearity  of  three  successive  control  points  does  guarantee  a 
continuously  varying  tangent.  However,  it  is  important  to  note  that  collinearity  of 
three  distinct  control  points  is  not  sufficient  to  guarantee  C'  continuity  as  explained  by 
Farin;  this  is  because  the  notion  of  C'  is  based  on  the  interplay  between  domain  and 
range  configurations.  Collinearity  of  three  points  is  purely  a  range  phenomenon. 
Without  additional  information  on  the  domain  of  the  curve  under  consideration,  we 
cannot  make  any  statements  concerning  differentiability.  However,  C‘  continuity  can 
exist  under  certain  conditions,  see  Farin  [Ref.  6:  p.  93]. 

The  distances  from  the  knot  point  to  adjacent  control  points  determines  the 
continuity  and  shape  of  the  curve  at  that  knot  point.  If  these  distances  were  equal,  the 
curve  would  have  continuity.  As  seen  in  Figure  6,  these  distances  are  usually  not 
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equal,  thus  the  curve  is  G'.  If  these  distances  become  zero,  or  the  control  points 
coincide  with  the  knot  point,  a  comer  will  be  formed.  This  will  be  shown  in  a  later 
example.  If  both  control  points  precede  or  follow  the  knot  point  a  cusp  will  occur. 
These  parameters  are  also  available  to  the  optimization  function.  Henceforth,  the 
control  points,  knot  points,  and  data  points  will  be  represented  as  in  Figure  6. 

D.  COMPUTING  NORMAL  DISTANCE  FROM  DATA  POINTS  TO  CURVE 
The  basic  problem  can  be  stated  thusly:  given  a  parametric  curve  F(t)  and  a 
point  Q,  both  in  the  plane,  find  the  point  on  the  curve  closest  to  Q.  In  other  words, 
find  the  parameter  value  t  such  that  the  distance  from  Q  to  F(t)  is  a  minimum.  Note 
that  the  line  segment  (whose  length  we  wish  to  minimize)  from  Q  to  F(t)  is 
perpendicular  to  the  tangent  vector  of  the  curve  at  F(t),  unless  the  closest  point  is  an 
endpoint  of  the  curve  segment.  The  equation  we  wish  to  solve  for  t  is 

[F(t)-Q]  •  F'(t)=0  (1) 

Curve  F(t)  is  a  cubic  Bdzier  polynomial  in  our  case; 

n 

F(c)  where 

are  the  control  points.  Curve  F'(t)  can  also  be  expressed  in  Bdzier  form; 

nt)  =  nT  . 

We  are  dealing  with  cubics,  so  polynomial  F(t)  is  of  degree  three,  and  F'(t)  is 
of  degree  two.  Polynomial  F(t)  -  Q  is  also  of  degree  three,  so  the  polynomial 
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described  by  equation  (1)  is  generally  of  degree  five.  Thus  the  problem  can  be 
restated  as  one  of  finding  the  roots  of  this  fifth-degree  polynomial.  There  is  no 
closed-form  solution  to  this  problem.  A  technique  developed  by  Schneider  [Ref.  19: 
pp.  607-611],  converts  the  equation  to  Bdzier  form,  and  then  uses  a  recursive 
algorithm  to  find  die  roots.  The  roots  are  evaluated  to  find  the  points  on  the  curve. 
By  comparing  the  distances  from  those  points  on  the  curve  to  the  arbitrary  point,  and 
also  considering  the  endpoints  of  the  curve,  the  desired  result  is  found  -  the  point  on 
the  curve  closest  to  the  arbitrary  point,  and  also  its  parameter  value. 

Plass  and  Stone  use  Newton-Raphson  iteration  to  find  the  root.  The  Newton- 
Raphson  iteration  for  solving  f(t)  =  0  is 


t 


fU) 

f{t) 


Referring  back  to  equation  (1),  each  iteration  will  reduce  t  by 

[fit)  -g]  'F'it)  ’t^ 

Because  Newton-Raphson  iteration  converges  quickly  and  is  fairly  inexpensive,  only  a 
few  steps  are  needed  to  find  a  close  approximation  to  the  root.  This  method  has 
proven  to  be  very  stable,  however  their  algorithm  is  not  suitable  for  fitting  smooth, 
continuous,  piecewise  functions  because  the  endpoints  and  tangents  for  each  cubic 
piece  cannot  easily  be  constrained. 

With  the  data  parameterized  by  normalized  accumulated  chord  length,  Marin 
and  Smith  [Ref.  20]  used  ODR  (Orthogonal  Distance  Regression)  splines  to 
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approximate  the  data.  Their  fit  is  better  as  compared  with  using  ordinary  least 
squares,  however  under  certain  circumstances,  there  will  not  exist  a  solution,  or  the 
solution  will  not  be  unique. 
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m.  IMPLEMENTATION 


A.  INITIAL  GUESSES 

The  algorithm  consists  of  several  MATLAB  programs  including  a  main  driver, 
”lsg2".  We  first  start  with  a  set  of  ordered  data  points,  input  by  the  user.  The  user 
can  either  select  the  knot  spacing  or  let  the  program  do  it.  Regardless,  the  initial  knot 
locations  are  a  subset  of  the  data  points.  The  program’s  knot  sequence  is  based  upon 
the  size  of  the  data  and  the  number  of  knot  points.  The  ultimate  goal  of  knot 
placement  is,  of  course,  to  use  as  few  knots  as  possible  to  get  a  "good"  representation 
of  the  shape  of  the  data. 

Next,  the  unit  tangent  vector  was  estimated  for  each  control  point.  This  was 
accomplished  by  using  a  chord  length  parametrization  to  fit  a  parametric  quadratic 
curve  to  five  data  points.  The  five  data  points  consisted  of  a  center  knot  point  and  its 
two  adjacent  data  points;  at  the  end  points  the  first  five,  or  last  five,  points  were 
used.  The  unit  tangent  vectors  were  approximated  by  the  unit  tangent  vectors  for 
those  parametric  quadratic  curves,  and  indicated  the  direction  of  the  adjacent  control 
points.  The  angles  of  the  unit  tangent  vectors  are  also  computed  as  they  are  the 
parameter  used  in  the  optimization  function. 

The  distance  from  the  knot  points  to  adjacent  control  points  were  calculated  by 
first  taking  one-third  of  the  distance  between  successive  knots.  The  unit  tangent 
vector  multiplied  by  that  distance  was  added  or  subtracted  from  each  knot  point  to 
obtain  the  locations  of  adjacent  control  points. 
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The  above  process  can  be  found  in  the  function  "iguess”.  The  data  points, 
knot  points,  and  control  points  are  then  plotted;  the  control  points  are  used  to 
calculate  the  points  of  the  approximating  cubic  Bdzier  curve.  Thus,  an  initial  guess  of 
the  shape  of  the  curve  is  obtained  and  plotted. 

B.  OPTIMIZATION  ROUTINE 

The  optimization  routine  used  was  "fmins”.  The  purpose  of  "fmins''  is  to 
minimize  a  function  of  several  variables.  For  example,  x=fmins(’fun’,xO)  returns  a 
vector  X  which  is  a  local  minimizer  of  fun(x)  near  the  starting  vector  xO;  ’fun’  is  a 
string  containing  the  name  of  the  objective  function  to  be  minimized.  In  our 
algorithm,  the  xO  array  consists  of  the  knot  points,  the  angles  of  the  unit  tangent 
vectors  and  the  distance  from  the  knot  points  to  adjacent  control  points.  These  are 
the  control  parameters. 

The  objective  function,  named  "objf2”  calculates  the  function  that  will  be 
minimized  by  the  optimization  routine  "fmins".  Along  with  the  x-array  which  was 
the  old  xO  array,  "objfZ"  also  calls  the  function  "ctpts"  that  computes  the  control 
points  and  "newk"  that  partitions  the  data  points  among  the  cubic  segments. 

Because  the  data  is  ordered,  it  is  necessary  for  the  closest  points  on  the 
approximating  curve  to  be  ordered  in  the  same  way.  Consequently,  it  is  necessary  to 
associate  the  data  points  with  particular  cubic  segments  to  avoid  computing  distances 
to  a  closer,  but  incorrect,  cubic  segment,  as  might  occur  if  the  data  makes  a  loop. 

The  function  "newk"  determines  which  data  points  partition  the  data  set  into  subsets 
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associated  with  the  various  cubic  segments.  These  data  points  will  be  called  dividing 
points.  So,  given  an  initial  knot  sequence,  the  search  for  the  dividing  data  points  for 
the  interior  knots  is  achieved  by  finding  the  smallest  distance  from  that  knot  point  to 
the  surrounding  data  points.  From  that  knot  point,  the  data  points  tested  are  those  up 
to  but  not  including  the  previous  and  subsequent  dividing  points.  For  the  first  and 
last  knots,  the  data  points  tested  are  those  up  to  the  next  and  previous  dividing  points 
respectively. 

During  the  optimization  routine,  the  knot  points  initially  coincide  with  the  data 
points,  and  the  subscripts  of  those  data  points  indicate  the  dividing  points  (which  are 
represented  by  the  global  variable  "dpkpc")  for  each  segment.  However,  as  the  knot 
points  move,  the  dividing  points  may  change.  With  each  new  iteration,  the  (possibly 
new)  dividing  points  must  be  obtained. 

When  the  function  "sod"  is  called,  the  sum  of  the  squares  of  the  distances 
from  the  data  points  to  the  nearest  point  on  the  cubic  segment  is  computed.  In 
addition,  the  square  of  the  distance  from  the  first  and  last  data  points  to  the  first  and 
last  knot  points,  respectively,  was  included  in  this  sum.  This  was  done  to  ensure  that 
the  curve  begins  near  the  first  knot  point  and  ends  near  the  last  knot  point. 
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IV.  CONCLUSIONS 


A.  EXAMPLES 

The  following  examples  illustrate  how  well  the  optimization  routine 
s^proximates  the  curve  given  the  initial  guess.  The  initial  guess  is  presented  first, 
followed  by  the  final  optimized  curve. 

A  major  limitation  of  the  initial  guess  curve  is  the  number  and  placement  of 
the  knot  points.  If  the  spacing  of  the  knot  points  is  too  close  together  the  initial  guess 
curve  will  make  a  straight  line  where  a  small  curve  should  appear.  Given  that  the 
curve  must  pass  through  the  knot  points,  if  the  data  points  or  knot  points  are  spaced 
far  apart,  certain  sections  of  the  curve  would  not  come  close  to  those  data  point,  s 
However,  that  distance  is  minimized  by  the  optimization  routine,  resulting  in  a  good 
approximating  final  curve,  as  will  be  seen  in  the  examples. 

Example  (1):  This  data  is  from  a  reacting  chemical  system,  taken  from  Marin 
and  Smith  [Ref.  20].  The  abscissa,  C^o.  of  each  data  point  is  an  input  concentration 
of  carbon  monoxide  in  a  catalytic  system.  The  ordinate,  R,  is  the  steady  state 
oxidation  rate  achieved  by  the  system.  The  initial  guess  curve  consists  of  one  interior 
knot,  and  is  not  a  "bad"  approximating  curve.  The  initial  knot  positions  are 
represented  by  "k"  for  the  initial  curve,  and  the  dividing  points  for  the  cubic  segments 
are  represented  by  "dpkpc",  which  in  some  sense  is  similar  to  new  knot  positions,  for 
the  final  curve.  Figures  7  and  8  shows  the  improvement  of  the  final  curve. 

Example  (2);  This  data  is  also  taken  from  Marin  and  Smith,  and  is  a  sample 
data  for  parametric  curve  in  R^,  with  seven  knots.  The  initial  guess  curve  is  not  very 
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Figure  7.  Initial  guess,  K=[l  13  23]. 
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good.  The  curve  loops  over  itself  and  thus  misses  data  points  by  great  distances,  see 
Figure  9.  But  the  power  of  the  optimization  function  reduced  this  distance  and 
removed  the  loop,  producing  an  excellent  approximating  final  curve.  Figure  10. 

Example  (3):  This  data  contained  eighty-two  data  points  with  eight  knots. 

The  data  represents  the  figure  of  the  letter  "H”,  Figures  11  and  12,  again  a  much 
improved  final  curve. 

Example  (4):  This  data  contained  thirty-six  data  points  for  the  curve  y  —  |  x  | . 
The  curve  contains  one  interior  knot  point  located  at  the  vertex.  The  data  forms  a 
straight  line,  but  the  initial  guess  curve  does  not,  because  of  the  distances  of  the  knot 
points,  see  Figure  13.  However,  in  Figure  14,  the  control  points  that  were  located 
near  x=±0.S  converged  to  the  knot  point  at  the  vertex,  forming  a  needed  comer. 
Thus,  the  final  curve  produced  was  an  exact  approximation,  to  within  the  convergence 
tolerance.  However,  note  that  the  "curve"  degenerates  to  two  straight-line  segments 
and  that  the  location  of  one  interior  control  point  is  arbitrary  on  the  line  segment. 

This  non-uniqueness  however,  did  not  affect  the  curve’s  convergence.  In  some 
examples,  such  as  (2)  and  (3),  the  maximum  iterations  allowed  by  "fmins"  was 
reached,  but  if  required,  the  continuation  of  iterations  is  incorporated  in  the  main 
driver  program  "lsg2". 

Figure  6,  which  contained  sixty-four  data  points  and  five  knots,  was  also  a 
very  good  final  curve  for  that  data;  the  initial  guess  curve  contained  no  loops  but 
nevertheless,  it  was  not  good. 
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Figure  11.  Initial  guess,  k=[l  9  20  38  49  59  68  82]. 
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lal  curve,  dpkpc=[l  9  16  37  49  67  82]. 
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Figure  14.  Final  curve,  dpkpc=[l  19  36]. 
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B.  CONCLUSIONS 


We  have  presented  a  method  for  approximating  ordered  data  by  using 
geometrically  continuous  piecewise  parametric  cubics.  After  choosing  the  knot  points 
which  were  a  subset  of  the  data  points,  we  used  local  quadratic  approximations  with 
chord  length  parametrization  to  estimate  the  unit  tangent  vectors.  In  turn,  the  unit 
tangent  vectors  are  used  to  calculate  the  position  of  the  interior  control  points  for  each 
cubic  segment.  Thus,  an  initial  guess  of  the  curve  was  approximated  using  the  above 
parameters.  The  optimization  routine  "fmins"  minimized  the  sum  of  the  squares  of 
the  distances  from  the  data  points  to  the  points  on  the  approximating  curve.  By  the 
examples,  this  optimization  routine  works  well,  however,  in  some  examples  such  as 
(3),  it  was  slow  to  converge  which  can  become  computationally  expensive. 

A  better  approach  to  this  problem  would  be  to  first  optimize  on  each  cubic 
segment,  which  can  be  achieved  by  fixing  the  knots  and  the  unit  tangent  vectors  for 
that  segment,  and  allowing  the  location  of  the  control  points  on  the  necessary  line  to 
be  free  parameters.  The  "sub”  optimization  routine  would  minimize  the  sum  of  the 
squares  of  the  distances  between  the  data  points  associated  with  that  cubic  segment 
and  the  curve.  Since  the  above  process  would  be  performed  for  each  cubic  segment, 
an  improved  initial  approximating  curve  would  result.  Choosing  a  better  initial  guess 
should  speed  up  convergence  to  the  final  approximation  curve. 

Another  approach  would  be  to  gear  the  optimization  process  to  this  problem 
instead  of  using  a  general  routine.  In  other  words,  an  algorithm  could  be  developed 
that  exploits  the  good  approximating  and  other  properties  of  piecewise  cubic 
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polynomials.  Exactly  how  this  can  be  done  remains  open,  but  the  process  is 
conceivable. 


A  large  volume  of  work  has  already  been  done  and  is  presently  continuing  in 
approximating  curves  to  fit  ordered  data  in  two  or  three  dimensions.  The  algorithms 
presented  here,  can  be  used  as  a  basis  for  further  work  by  geometric  modelling 
researchers  and  graphics  programmers  in  the  field  of  computer  aided  geometric 
design. 
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APPENDIX:  MAIN  HtOGRAMS 


Pn^iram  Isg2 


%  Program  lsg2 

%  This  program  is  the  main  driver  for  Least  Squares  Approximation  By  G1 
%  Piecewise  Parametric  Cubics.  It  takes  data  points, Q;  and  uses  the 
%  the  function  'iguess'  to  get  a  plot  of  the  initial  guess  curve.  It  also 
%  uses  the  array  returned  by  the  function  'iguess'  for  input  to  the 
%  optimization  function  'FMINS' .  Furthermore,  it  enables  the  user  to  continue 
%  with  the  iterations,  if  the  maximum  iterations  are  reached  by  the 
%  optimization  routine. 

global  dpkpc 

disp( 'Choose  one  of  the  following:  dpts,  dptsl,  ...,  dptsN' ) 

Q  =  input ('  '); 

[r,m]  =  size(Q) ; 

%  Error  checking. 

if  r  ~=  2 

disp('Set  of  points  must  be  2  x  n  matrix,  try  again.'),  pause (2) 
lsg2 

end 

%  Obtains  the  xO  array  for  input  to  the  optimization  function  'FMINS'. 
xO  =  iguess (Q ) ; 
pause (5) 

options  =  [0 , 1 . e-3 , 1 . e-3 ] ; 
tb  =  clock; 

%  Calls  the  function  'FMINS'. 

X  =  fmins ( 'objf2 ' ,  xO,  options,  [],Q); 
et  =  etime (clock, tb) ; 

%  Enables  the  user  to  continue  with  the  iterations,  start  anew  or  end. 

disp("rype  "1"  to  continue  iterations;  "2*  to  start  anew,  "3''  to  end') 
b  =  input ( '  ' ) ; 


49 


if  b  ==  1 

disp('  How  many  more  iterations  would  you  like? 
o  =  input ( '  ' )  ; 

options (14)  =  o; 
xO  =  x; 

X  =  fmins { 'obj f2 ' ,  xO,  options,  [),  Q)  ; 

elseif  b  ==  2 
lsg2; 


elseif  b  ==  3  I  b  -=  1  I  b  -=  2  ; 

disp{'  This  progreim  will  now  end!  ') 

end 
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function  G  =  iguess(Q) 


%  function  G  =  iguess(Q) 

%  This  function  takes  data  points,  Q;  and  with  a  subset  of  the  data  points 
%  computes  knot  points,?;  the  position  of  the  knot  points, k;  and  the  distance 
%  between  the  successive  knot  points, dt.  It  calls  the  functions  'unitv'  to 
%  compute  the  unit  tangent  vectors  for  each  knot  point,  'ctpts'  to  compute 
%  the  control  points,  and  'plotC'  to  plot  the  initial  guess  curve.  The  function 
%  returns  the  array  that  contains  the  initial  guess  curve's  parameters.  This 
%  array  will  be  used  by  the  main  driver  progreun  'lsg2'  for  input  to  the 
%  optimization  function  'FMINS'. 


global  dpkpc 
[r,m]  =  size{Q) ; 

disp('  Give  the  number  of  knotpoints') 
n  =  input ( '  ' ) ; 

disp('Type  ‘‘I"  for  default  knot  position  or  “Z"  to  input  your  own.') 
h  =  input ( '  ' )  ; 

if  h  ==  1 

k  =  round( ( (m-1) / (n-1) ) * [0 ;n-l]  +  on6s(l,n)); 
k; 

elseif  h  ==  2 

disp( 'Input  the  initial  knot  point  position,  example;  [148  ...n]') 
k  =  input ( '  ' )  ; 

elseif  h  ~=  1  I  h  ~=  2 

disp('  You  can  only  choose  *1''  or  *2".  You  must  start  over .'),  pause  (2 ) 
iguess 

end 

dpkpc  =  k; 

%  Gets  the  knotpoints 
P  =  [P  Q( k) ] ; 

%  Computes  the  distance  between  successive  points 
[s,t]  =  size(P) ; 
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d  *  P(:,l:t-1)  -  P(:,2:t); 

d  a  (1/3)  *  sqrt  (suin(d.  *d) )  ; 
dt  =  [d;d]; 

%  Calls  the  function  "unitv*  to  compute  the  unit  tangent  vectors, 
u  =  unitv(Q, k) ; 

ang  =  atan2 (  u ( 2 , : )  , u (1 , : )  ) ; 

%  Calls  the  function  "ctpts*  to  compute  the  control  points. 

C  =  ctpts (P, ang, dt) ; 

%  Calls  the  function  "plotC"  to  plot  the  initial  guess  curve, 
plot  =  plotC(C,Q, P) ; 

%  Sets  up  the  array  that  contains  the  initial  guess  curve's  pareuneters 
%  that  are  available  to  the  optimization  routine. 

G=  [P(l,:)  P(2,:)  ang  dt(l,:)  dt{2.:)]; 


function  os  =  objf2(x,Q) 


%  function  os  =  objf2(x,Q) 

%  This  is  the  objective  function  that  will  be  minimized  by  the  optimization 
%  function  'fmins'.  The  input  arguments  are  the  vector  x  that  minimizes 
%  function  objf2(x),  and  the  data  points  Q.  The  output  is  the  sum  from 
%  the  function  'sod'  plus  the  distance  squared  from  the  first  and  last  data 
%  points  to  the  first  and  last  knot  points. 

%  Begin  function 
global  dpkpc 
m  =  length (x) ; 
n  =  round (m/ 5) ; 

[r, s]  =  size(Q) ; 


%  Picks  out  the  knot  points. 

P(l, :)  =  X(l:n) ; 

P(2. : )  =  x{n+l:2*n) ; 

%  Picks  out  the  angles  of  the  unit  tangent  vectors, 
ang  =  x(2*n+l ; 3*n) ; 

%  Picks  out  the  distances. 
dt(l,:)  =  x(3*n+l : 4*n-l ) ; 
dt(2, : )  =  x(4*n:m) ; 

%  Calls  function  that  compute  the  control  points. 

C  =  ctpts (P, ang, dt) ; 

%  Calls  function  that  computes  the  new  dividing  point  positions, 
dpkpc  =  newk(Q,P); 

%  Calls  the  function  that  computes  the  sums  of  the  square  of  the  distances 
%  from  the  data  points  the  nearest  point  on  the  cubic  segment.  The  distance 
%  from  the  first  and  last  data  points  to  the  first  and  last  optimization 
%  points  are  also  added  to  this  sum. 

fp  =  P(:,l)  -  Q{:,1); 

Ip  =  P{ : ,n)  -  Q( : , s) ; 


os  =  sod(C, Q, dpkpc)  +  fp'*fp  +  lp'*lp; 


disp(  [  os  ] ) 
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function  uv  =  unitv(Q,k) 


%  function  uv  =  unitv(Q,k) 

%  This  function  takes  data  points  'Q'  from  a  separate  file,  and  the  position  of 
%  knotpoints  'k'as  input  variables.  It  uses  chord  lengths  parameterization  to 
%  fit  a  parametric  quadratic  curve  to  five  data  points.  The  unit  tangent 
%  vectors  are  approximated  by  the  unit  tangent  vectors  for  these  quadratic 
%  functions .The  output  is  the  set  of  unit  tangent  vectors  in  the  direction  of 
%  the  knot  points . 

%  Begin  function. 

[r.m]  =  size(Q) ; 

n  =  length ( k) ; 

for  j  =  l:n 

if  j  ==  1,  k(j)  =  1;  kt  =  1; 
elseif  j  ==  n,  k(j)  =  m-4;  kt  =  5; 
else  k(j)  =  k(j)-2;  kt  =  3; 
end 

%  Extracting  the  knot  point  and  four  adjacent  points 
X  =  Q{l,k( j) :k( j)+4) ' ; 
y  =  Q(2,k(j) :k(j)+4) 

%  Calculating  the  chord  lengths 
xd  =  dif f (x) ; 
yd  =  diff (y) ; 

d  =  sqrt (  xd.*xd  +  yd.*yd); 
t(l)  =  0;  t(2)  =  d(l); 

t{3)  =  t(2)  +  d(2); 

t(4)  =  t(3)  +  d{3) ; 

t(5)  =  t(4)  +  d(4); 

%  Approximating  coefficients  of  quadratic 

c  =  [ones(5,l)  t'  (t.*t)']  \  [x  y] ; 

u  =  c(2, :)  +  2*c(3, :)*t(kt) ; 

u  =  u  /  norm{u) ; 

uv ( ; , j )  =  u ' ; 


end 
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function  nk  =  newk(Q,P) 


%  function  nk  =  newk(Q,P) 

%  This  function  takes  data  points,  Q;  and  knot  points,  P;  as  input.  The 
%  function  finds  the  closest  data  point  of  that  cubic  segment  that  is 
%  associated  with  that  knot  point,  and  returns  a  new  k-array,  nk. 

%  dpkpc  is  a  global  variable  that  is  initially  equal  to  the  old  k-array,  k. 

%  Begin  function. 

global  dpkpc 

[r,ro]  =  size(Q) ; 

[s, n]  =  size (P) ; 

nk(l)  =  1;  nk(n)  =  m; 


for  i  =  2;n-l 

js  =  dpkpc (i-1);  je  =  dpkpc (i+1);  jm  =  dpkpc (i); 
z  =  je-js+1;  mm=jm-js+l; 

R  =  Q(;,js:je)  -  P(:,i)  *  ones(l,z); 

for  j  j  =  1 :  z 

D(jj)  =  R(:,jj)'  *  R(:,jj); 

end 

sd  =  sign(  D(min)  -  D(mm+1)  )  ; 
while  D(mm)  -  D(mm+sd)  >  0 
if  mm  ==  2  &  sd  <  0 
break,  end 

if  mm  ==  m-1  &  sd  >  0 
break,  end 

mm  =  mm  +  sd; 
end 

nk(i)  =mra+  js  -  1; 

end 

dpkpc  =  nk; 


55 


function  sumofdist  =  sod(C,Q,dpkpc) 

%  function  sumofdist  =  sod(C, Q, dpkpc) 

%  This  function  takes  control  points. C;  data  points, Q;  and  the  position 
%  of  the  data  point  closest  to  that  knot  point  within  that  segment.  The 
%  function  returns  the  sum  of  the  distance  squared  from  the  data  points 
%  to  the  nearest  point  on  the  cubic  segment . 

%  Begin  function 

n  =  length (C) ; 

[r, s]  =  size(Q) ; 

y  =  dpkpc; 
t  =  length (dpkpc) ; 

%  Initialize  counter  and  sum. 

cntr  =  0 ; 
sum  =  0 ; 

for  i  =  l:3:n-3 

cntr  =  cntr  +  1; 

for  j  =  y(cntr) :y(cntr+l) 

xy  =  NearestPoint (  C(:,i:i+3)',  Q(:,j)');  %  See  note, 

d  =  (  Q(:,j)'  -  xy  ); 
sum  =  sum  +  d  *  d ' ; 
if  j  ==  y(cntr)  &  i>l 
d2  =  d*d'; 
ds2  =  ds*ds'; 
dm  =  max(d2,ds2); 
sum  =  sum  -  dm; 

end 

end 

ds  =  d; 

end 


sumofdist  =  sum; 

%  Note:  NearestPoint;  obtained  from  'Solving  the  Nearest  Point-on-Curve 
%  Problem'  and  'A  Bezier  Curve  Root-Finder'  developed  by  Philip  J. 

%  Schneider  in  “Graphics  Gems“,  Academic  Press,  1990. 
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functioo  C  s  ctpts(P,anig,dt) 


%  function  C  =  ctpts ( P, ang, dt ) 

%  This  function  takes  knot  points, P;angle  of  the  tangent  vector,  ang  ;  and 
%  the  distance  between  successive  knot  points, dt;  as  input.lt  then 
%  computes  the  unit  tangent  vector  and  the  control  points. 

%  Begin  function. 

n  =  length (P) ; 


%  Converts  the  angle  to  its  x  and  y  components, and  computes  the  control  points. 


for  k  =  2:n-l 

u  =  [cos (ang (k))  ;  sin (ang (k) ) ] ; 

T  =  [T  P(:,k)-u*dt(2,k-l)  P(:,k)  P ( : , k) +u*dt ( 1 , k) 1 ; 

end 

ul  =  [cos(ang(l) ) ;  sin(ang(l) ) ] ;  un  =  [cos(ang(n))  ;  sin (ang (n) ) ) ; 
C=  [P(:,l)  P(:,l)+ul*dt(l,l)  T  P ( : , n) -un*dt ( 2 , n-1 )  P(:,n)]; 


57 


function  plotC  »  pltC(C,Q,F) 


%  function  plotC  =  pltC(C,Q,P) 

%  This  function  takes  as  input:  control  points, C;  data  points, Q;  and  knot 
%  points.  The  control  points  are  used  to  calculate  the  points  of  the 
%  approximating  cubic  Bezier  curves.  The  control,  data,  and  knot  points 
%  are  then  plotted  along  with  the  curve. 

%  Begin  function. 

[s, t]  =  size(C) ; 

X  =  [0:. 025:1]; 

[a,b]  =  size(x) ; 

W  =  [  ]  ; 

for  j  =  l:3:t-3 

Y  =  zeros (2, b) ; 

M  =  (berny(3, 0,x) '  berny (3, l,x) '  berny(3, 2,x) '  berny ( 3 , 3 , x) ' ] ; 

Y  =  Y  +  C(:,j:j+3)  *  M'  ; 

W  =  [W  Y]; 

end 

plot(  W(l,:)  ,  W(2,:)  ),  hold 
plot(  C(l, :)  ,  C(2, :)  ) 

plot(  Q(l,  :)  ,  Q(2,  :)  ,  '  +  ') 

plot(  P(l, :)  ,  P(2, :)  ,  'X') 

plot(  C(l, :)  ,  C(2, :)  ,  'o') 
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function  pop  =  poplt(x,Q) 


%  function  pop  =  poplt(x,Q) 

%  This  function  picks  out  the  knot  points,  angles,  and  distances  from 
%  the  optimization  function  x  =  'fmins'.  It  then  calls  the  function  that 
%  computes  the  control  points.  Using  the  given  data  points;  Q,  a  plot  of 
%  the  piecewise  cubic  curve  is  returned. 

%  Begin  function. 

m  =  length (x) ; 
n  =  round (m/ 5); 

%  Picks  out  knot  points,?;  angles, ang;  and  distances,  dt . 

P(l,:)  =  X{l:n);  P{2,:)  =  X(n+l:2*n); 

ang  =  x(2*n+l : 3*n) ; 

dt(l,:)  =  x(3*n+l;4*n-l) ;  dt(2,:)  =  x(4*n:m) ; 

%  Calls  the  function  that  computes  the  control  points. 

C  a  ctpts (P, ang, dt) ; 

%  Calls  the  function  that  plots  the  piecewise  cubic  curves. 
plotC(C,Q, P) 
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function  val  =  berny(n,  i,  x) 


%  function  val  =  berny(n,  i,  x) 

%  This  function  is  a  non-recursive  formula  for  Bernstein  Polynomials  which 
%  form  a  basis  for  Bezier  curves.  The  inputs  are  the  degree  of  the  poly- 
%  nomial,  n;  the  particular  curve  that  is  assigned  a  value  of  zero  up  to 
%  and  including  the  degree,  i;  and  the  points  between  [0,1]  to  be  evaluated,  x. 
%  The  output  is  the  coordinates  of  points  on  the  curve. 

%  Begin  function. 

ni  =  [1  3  3  1] ;  %  See  note. 

m  =  size (x) ; 


if  n  <  i 

val  =  zeros (m); 
eiseif  i  <  0 

val  =  zeros (m); 
eiseif  ( (n  ==  0)  &  (i  ==  0)) 
val  =  1 ; 
else 

val  =  ni(i+l)  *  (x.^'i)  .*  {(ones(m) 

end 


X)  .''(n-i)); 


%  Note;  This  function  will  only  work  for  cubics  (which  are  used  throughout 
%  the  supporting  programs) .  For  that  reason,  it  is  more  efficient  to  use 
%  the  coefficients  for  a  third  degree  polynomial  rather  than  one  for  a 
%  general  n-degree  polynomial . 
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